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Abstract 

Ah initio density functional theory has been used to analyze flexural modes, elastic constants, 
and atomic corrugations on single and bi-layer graphene. Frequencies of flexural modes are sensitive 
to compressive stress; its variation under stress can be related to the anomalous thermal expan- 
sion via a simple model based in classical Elasticity Theory.l' Under compression, flexural modes 
are responsible for a long wavelength rippling with a large amplitude and a marked anharmonic 
behavior. This is compared with corrugations created by thermal fluctuations and the adsorption 
of a light impurity (hydrogen). Typical values for the later are in the sub- Angstrom regime, while 
maximum corrugations associated to bending modes quickly increase up to a few Angstroms under 
a compressive stress, due to the intrinsic instability of flexural modes 
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I. INTRODUCTION. 



Physical properties of single and multilayered graphene are remarkable and open interest- 
ing avenues for applications in the futureP However, to realize its potential it is important 
to fully understand the elastic properties of graphene layers^ the effect of external stresses 
on electronic propertieP'^, or conversely the effect of doping on mechanical properties,'^ and 
the appearance of ripplings,'^ corrugations, and other defects expected under actual con- 
ditions of cleanness, temperature, and growth. The understanding of any departure from 
the perfect 2D configuration (e.g., rippling and corrugation) of very thin graphene layers, 
either in free suspension or deposited on a template substrate, is central for manufacturing 
processes.'^'^ At T= K, in the absence of defects, the carbon bond on the graphene layer 
is well understood in terms of the formation of three in-plane localized strong sp^ bonds, 
and a fourth delocalized, out-of-plane, vr-like bond.'^'^ The optimum geometrical configu- 
ration is achieved by a honeycomb lattice formed by two equivalent sublattices displaying 
P6/mmm symmetry. The corresponding electronic structure shows bands dispersing linearly 
around the Fermi energy that are responsible for the fast and efficient transport of carri- 
ers. Both experimentally and theoretically,'^'^^'^ it is shown that this kind of arrangement 
results in a material with the largest in-plane elastic constants known yet. Therefore the 
2D perfect flat layer makes the most stable configuration since deviations from a common 
plane requires a significant amount of energy. Any departure from such a scenario affects 
greatly the atomic scale properties of the layer and must be understood in order to effi- 
ciently exploit graphene's properties. Different reasons, however, might be invoked for a 
two-dimensional graphene layer to adopt a certain corrugation at different scales. First, at a 
non-zero temperature there is a thermodynamic argument that implies the impossibility for 
a perfect 2D layer to exist in SD.EEHIIHIZI Second, defects like adsorbed impurities, vacancies, 
etc, create local corrugations at the atomic scal^ that propagate via the elastic properties 
of the lattice originating long-range correlations. Finally, external applied stresses related 
to conditions on the boundary make graphene to bend and to corrugate; an interesting 
point to study since the growth of graphene layers on different supporting substrates im- 
plies mismatches that introduce all kind of stresses that have been observed to originate a 
highly complex and corrugated landscape. '^^'^ Research on these issues require atomically 
resolved techniques; theoretical simulations performed with ab-initio Density Functional 
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Theory (DFTp^^ give an accurate picture of the interatomic interactions and can provide 
valuable information about geometry, vibrations, electronic effects, etc, down to the atomic 
level. Indeed, ab-initio DFT has been successfully applied for many years now to describe 
the different allotropes of carbon, like graphite, diamond, fullerenes, organic molecules, nan- 
otubes, etc (e.g. see refsJ ^^*^^ *^. Here, we focus on the role of external stresses on the 
flexural modes ultimately responsible for the intrinsic thermodynamic instability of these 
layers. The frequencies of these modes are very sensitive to external compressive stresses and 
can be related, in turn, to the Griineisen coefficients and the thermal expansion coefficient 
of graphene.'^ The purpose of this paper is to use state-of-the-art ab-initio DFT to compare 
atomic corrugations on graphene layers related to three different origins: (i) the presence 
of small adsorbed impurities (hydrogen), (ii) thermal vibrational amplitudes (T=300 K), 
and (iii) the long wavelength corrugation characteristic of flexural modes. We analyze the 
consequences of an isotropic compressive stress on these corrugations, the elastic constants 
and the phonons. 



II. AB-INITIO DENSITY FUNCTIONAL ATOMISTIC MODEL. 

Since we are interested in elastic and vibrational properties we use accurate norm- 
conserving pseudopotentials to describe carbon valence band electrons (2s^ 2p'^)P^^ Elec- 
tronic wavefunctions have been expanded in a plane-wave basis set up to an energy cutoff 
of 750 eV, and a 12 x 12 x 2 Monkhorst-Pack mesh has been used to sample wavefunctions 
inside the first Brillouin zone. Electronic bands were obtained using a smearing width of 
rj = 0.1 eV. The choice of the exchange and correlation (XC) potential opens up various 
possibilities. We favor the local density approximation (LDA)PSI because of its simplicity, 
excellent structural results for a wide range of carbon allotropes relevant for this work, and 
in particular because, even if somehow fortuitous, it predicts a reasonable value for the layer 
separation in graphite. LDA presents a solid record reproducing many key features for a 
wide variety of carbonaceous materials, and it has been a standard choice in the literature 
for many authors. Unfortunately, the use of LDA cannot be justified "a priori" on physi- 
cal grounds by merely referring to the behavior of the next term on a gradients expansion 
(GGA); different flavors of GGA yield worse agreement with experiments (e.g. for the dis- 
tance between planes in graphite) showing that the converge on such an expansion is not a 
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trivial one. Therefore, we are forced to justify the use of LDA on the "a posteriori" reason- 
able structural results. Its simpler formulation should also be taken as an advantage, since 
makes a better defined and more unique formulation. Finally, it is known that LDA may 
not be too accurate for absolute values of total energies, but it affords a good description of 
relative values (e.g., barriers), and it yields the correct dependence with stress after a single 
common offset correction.l^Zl A periodic supercell with a vacuum gap in the perpendicular 
direction (30 A) and several n x m in-plane periodicities (?T,,m = 1 — 8) depending on the 
requirements of each model have been used. Geometrical parameters have been optimized 
so residual maximum forces and stresses have been kept below F^ax < 0.0001 eV/A, and 
Smax < 0.0002 GPa for the calculation of phonons and elastic constants, yielding a lattice 
parameter for the rhombohedral unit cell on a single graphene layer of a = 6 = 2.447 A 
(Fig. [l]), and c = 6.65 A for graphite (Bernal stacking). Calculations were performed with 
the CASTEP code, as implemented in Materials Studio.'^^'^ 

The elastic tensor has been obtained by applying a series of strains on the 1x1 unit cell 
(maximum amplitude 0.003 A), and by computing the associated stress tensors Phonons 
have been computed using linear response theory.'^ As a global benchmark to the elements 
included in our calculations (pseudopotential, exchange and correlation model, different 
convergence thresholds, etc), we compare in Table I our theoretical results for the five in- 
dependent components of the elastic tensor for graphite to experimental values measured 
on pyrolitic graphite by ultrasonic resonance,'^ and on single crystals by inelastic x-ray 
scattering.'331 A comparison to a recent independent theoretical determination has also been 
included.®! The agreement is quite reasonable, except perhaps for cis that cannot be de- 
scribed properly with a local choice for exchange and correlation. We notice that the 
basal in-plane theoretical Poisson ratio, u^y = Vy^ ~ ^ ~ 0.19, and the Young modu- 
lus, Yx = Yy = 1 (TPa), are in good agreement with experimental values, 0.165 and 1.1 TPa 
respectively!^ as expected these are very similar for a graphene monolayer (0.18 for theory 
to be compared with a quoted experimental value ofQ.l'^. 

Values for the 2D elastic constants of a single graphene layer (I),'^the bilayer with the 
usual AB Bernal stacking, and a bilayer with direct AA stacking (distance between layers 3.5 
A) display little variation in the main elements of the elastic tensor (Table II, for reference 2D 
tensions for graphite are quoted). The main departure has been found for the AA bilayer in 
the metastable configuration (distance between layers 1.5 A)Slthe reduction of Cn is related 
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to the gradual transformation of the in-plane C-C bond order from sp^ towards sp^, since the 
2D unit cell needs to be expanded to make the structure metastable under its own pressure. 
The increase in cge implies a larger tendency against shear due to substitution of the weak 
van der Waals-like interaction between both layers by stronger sp^-like chemical bonds. 

For the two-dimensional graphene layer, due to its 6-fold symmetry, we only need to 
determine two independent values in the elastic tensor; these can be related to the two 
Lame coefficients needed to characterize an isotropic system; A = Cyi and /i = Cge- An 
independent determination of Cn, that should be equal to C12 + 2c66, holds fairly well in all 
the cases. Therefore, simple analytical models based in the theory of elasticity that only 
need values for these two parameters make a good approximation and can be thoroughly 
analyzed.l^^ Under a moderate compressive stress (e.g. e = —0.04) the main change happens 
in the value of C12 that is reduced by ~ 2. 

Fig. [2] gives the phonons for graphene along two main directions in the irreducible part 
of the Brillouin Zone: G-M and G-K (Fig. [T]). Since there are two atoms in the unit cell, 
we find three acoustic (ZA, TA and LA) and three optical modes (ZO, TO and LO). The 
optimized structure shows a computed vibrational spectra that agree well with previous 
calculations in the literature,!^ as can be seen by comparing frequencies for the modes at 
selected high-symmetry points (Table III). For e = (no external stress, black circles), all 
the frequencies are positive and behave as expected; it is interesting to observe how the 
lowest acoustic mode (ZA) disperses quadratically, a characteristic feature for a 2D system 
(e.g., see Eq. (3) in re f.lil). 

A. Instability of flexural modes under compression. 

The intrinsic instability of the layer against compressive strain is demonstrated by the 
appearance of negative frequencies in a small region around the G point; the ZA phonon 
softens and goes to negative values.'^^'^ It is worth noticing that because the extreme asym- 
metry between compression and tension for strain induced instabilities on graphene layers 
we have only considered compressive strains.!^ Other failure mechanisms for graphene un- 
der tension have been studied elsewhere.!^ For some critical wavevector, q = Qc, the mode 
becomes again positive, and the value of Qc increases with the amount of stress. The value 
of Qc defines a characteristic wavelength, Ac = with an amplitude that depends directly 
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on the strain. This lowering of the frequency of the ZA modes under compression is consis- 
tent with the negative Griineisen coefficients reported in the literature J ^ * ^^ * ^^ * "^^ A fit to the 
results in Fig. [ijgives ix'l(THz^) ^ —25 |q|^ + 112 |q|^ for e = 0.02 and |q| in A^^ (valid near 
the center of the Brillouin Zone). For a -2% strain the corresponding 2D stress tensor is 
o" = CTxx = CTyy = 0.6 cV/A^, Qc ~ 0.35 A~^, and Ac ~ 12 A. Other relevant feature under 
stress observed in Fig. |2]is the hardening of LO and TO modes by ~ 20% for e = —0.04. 
These modes correspond to displacements inside the plane of the layer, play an important 
role in the dissipation of energy after the absorption of electromagnetic radiation ,S3E3 ^nd 
are strictly degenerate at the G point. 

Note, finally, that the softening of the ZA and ZO modes under strain, shown in Fig. [2] 
will enhance the intrinsic spin-orbit coupling of graphene.l^SI 

III. CORRUGATIONS UNDER STRESS. 
A. Flexural modes 

We address the question of deformations related to the intrinsic instability brought by 
the fiexural AZ mode. Under compressive strain, the negative region of frequencies near 
the G point indicate the existence of an energetically favorable deformation for the system. 
Therefore, we essay a cosine-like perturbation with a wavelength corresponding to qc on a 
rectangular 8x2 supercell under a range of strains from -2% to -5%. These strains are 
large enough for the wavelength determined by the 8x2 supercell. It has been previously 
shown by Zhang and Liii^ that wavelengths of periodic undulations scale inversely with the 
square root of strain, setting a minimum value for a given wavelength that it is < 1% for 
the 8x2, therefore compatible with the values we use here. For e =-3%, the soft mode 
displayed in Fig. |3] allows the system to gain —26 meV with a maximum amplitude of cm = 
max(zj) — min(zj) = 2.03 A, and an averaged corrugation of c = ^i=i n I ^^^i H 0-66 A. 
Larger stresses drives the non-linear behavior of the system, for e =-5% we find an energy 
gain of —3059 meV and cm = 2.81, c = 0.92 A. It is interesting to notice that the density of 
states at the Fermi level (e.g., the Tersoff-Haman STM image) for such a long-range rippling 
of graphene is symmetric with respect to the two sublattices, as can be seen in Fig. |3| where 
the density of states has been drawn on each C atom for 1^ = 1 V. Corrugations of this 
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kind are expected to evolve smoothly into long wavelength deformations in larger unit cells, 
unlike the case for thermal or defect induced corrugations we study below. 

B. Thermal fluctuations and light impurity adsorption. 

To quantitatively compare corrugations on graphene layers originated from other effects, 
we investigate the effect of thermal fluctuations (T 7^ K). The relationship between stresses 
originated in free finite edges and thermal fluctuations has already been reported in the 
literature.^®. In our simulations, based in periodic models, stresses due to free edges do not 
naturally occur, but it is interesting to notice how our uniform compressive strain of —0.05 
yields results comparable to those simulations. To introduce a finite temperature in our 
simulations we use ab-initio molecular dynamics to follow the trajectories of carbon atoms 
in a 4 X 4 unit cell in equilibrium with a Nose-Hoover thermostat at T = 300 K. Starting 
out from an equilibrium distribution the graphene layer is allowed to evolve in the canonical 
ensemble for tq = 1 ps using a At = 0.5 fs step to integrate the equations of motion. As 
before, the same two statistical indicators, c and cm-, are computed averaging over all the 
configurations reached inside the time interval tq (the shape of the layer at a time chosen 
at random, t = 260 fs, is shown in the upper panel of Fig. |4]). Table III shows values for 
three strains on the unit cell: e = 0, —0.03, and —0.05. The last two strains correspond to 
stresses of axx = CTyy = —0.94 eV/A^, and —1.69 eV/A^ respectively. The case e = allows 
a double check by comparing with the experimental Debye- Waller factor determined from 
low-energy electron diffraction on graphite (0001):'^'' 0.053 A to be compared with a value 
of 0.069 A as derived from the simulation. Even for the small strain values considered here, 
it is clear the non-linear growth of corrugations vs the external strain. 

C. Adsorption of a light impurity. 

A further source for corrugations at the atomic level is the presence of defects. For the 
sake of simplicity, we consider the adsorption of a light element (hydrogen atom), which is 
known to originate a static distortion in the lattice that depends on the coverage, 6, and 
extends as a long-range elastic perturbation.'23 For 6 = j^, a value where the interaction 
between periodic images of the adsorbate is already low and can be representative of the 
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behavior for the isolated impurity, both the averaged and the maximum corrugation take 
values fairly similar to the ones obtained by heating the layer to T = 300 K, except perhaps 
for the value of cm, that systematically takes lower values for the thermal case (Table 
III). The corrugations studied here lead to partial sp^ hybridization of the carbon orbitals, 
enhancing the spin-orbit coupling.'^ 



IV. CONCLUSIONS. 

We find that averaged atomic corrugations of ~ 0.05 A, and maximum corrugation values 
of ^ 0.5 A are easily obtained under realistic conditions of cleanness and temperature on 
graphene layers. Corrugations due to thermal vibrations of atoms at room temperature 
make a similar effect to the adsorption of a light atom (H) with a moderately small coverage 
{9 = 1/16). An external compressive stress increases corrugations in a non-linear way. Due 
to the flexural mode, the graphene layer becomes intrinsically unstable and shows a value 
for the maximum corrugation under compressive stress three or four times larger than values 
attained at T < 300 K or because the adsorption of impurities. Upon a 4% contraction, the 
LO and TO modes increase their frequencies by about 20 %; while the main elastic constant 
affected is Ci2, that is approximately halved.'^S'^ On the other hand, the frequencies of the 
ZA and ZO modes are reduced. These modes induce an effective sp^ coordination in the 
lattice, and their softening under strain leads to an enhancement of the intrinsic spin-orbit 
couplin^3SI^ opening a way to observe topological insulator features in graphene.!^ Finally, 
these simulations have been useful to obtain the elastic constants needed to construct an 
analytical model of free standing membranes based in elasticity theory.'^ Such a simple 
analytical model yields negative Griineisen parameters related to the bending modes, and a 
simple explanation for the negative thermal expansion coefficient of graphene. 
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TABLE I. Comparison of ab-initio elastic constants (GPa) computed for graphite with experimental 
values determined on single crystals by inelastic x-ray scattering'^^, on highly oriented pyrolitic 
graphit^^, and independent theoretical calculationd^D^ 
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TABLE II. Elastic constants (eV/Ang^) normalized to the carbon mass in the 1x1 graphene unit 
cell. Results are given for (I) a single graphene layer, (AB) bilayer with Bernal stacking, (AA) 
bilayer with direct stacking, (AA') meta-estable bilayer at short distanced, (BULK) graphite, 
and experimental values for graphite (EXF^). DFT calculations have been performed using LDA 
and norm-conserving pseudopotentials as explained in the text. No attempt to compute elastic 
constants involving strains in the z-direction has been made for super-cells with a vacuum separator. 
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TABLE III. Comparison of selected phonon frequencies (THz) at high-symmetry points in the 
Brillouin zone. 
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e = 0.00 


e = -0.03 


e = -0.05 


CASE nxn 


C CM 


C CM 


C CM 


T=300 K 4x4 


0.069 0.572 


0.153 1.280 
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H 4x4 
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0.310 1.580 


PHONON ZA 8x2 


0.000 0.000 


0.662 2.034 


0.917 2.807 



TABLE IV. Average, c (A), and maximum, cm (A), corrugations in a graphene layer subject to an 
external thermal bath (T = 300 K), the adsorption of a light atom (H), or related to the intrinsic 
instability due to its 2D character (mode AZ in^. 






FIG. 1. (Color online) Graphene rhombohedral 1x1 unit cell (7 = 60), and corresponding 
Brillouin Zone (special points are labelled according to the phonon calculation in Fig. [2]). 
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FIG. 2. (Color online) Phonons for the free-standing graphene layer (THz). Path: M(0, 2,0) — >■ 
G(0,0,0) K{1, 1,0). Black, circles: e = 0. Red, stars: e = -0.02. Blue, squares: e = -0.04. 
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FIG. 3. (Color online) Long-wavelength deformation related to the AZ mode computed on a 
8x2 rectangular unit cell (e = —0.03). Average and maximum corrugations are, c = 0.66, and 
Cm = 2.03 A, respectively. The density of states at Ep + 1 V over the C atoms is shown (STM 
image in the Tersoff-Haman approximation). 
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FIG. 4. (Color online) Upper pannel: Snapshot configuration for a graphene layer subject to 
vibrations at T = 300 K (an arbitrary time, t = 260 fs has been chosen). Lower pannel: equilibrium 
positions on the graphene layer upon hydrogen adsorption. 
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